##########################################################################################################
#R Script to replicate the results in:
#Model Selection in Observational Media Effects Research: A Systematic Review and Validation of Effects
#http://dx.doi.org/10.1080/00323187.2017.1411766
#For questions please contact Martijn Schoonvelde, mschoonvelde@gmail.com
##########################################################################################################

rm(list = ls())

library(memisc)
library(foreign)
library("arm")
library(ggplot2)
library("MatchIt")
library(plyr)
library(dplyr)
library("Zelig")
library("ZeligChoice")
library(xtable)
library(stargazer)
library("reshape")
library("apsrtable")
library(Hmisc)
library(modeest)
library(ordinal)


#remove scientific notation
options(scipen = 999)
options(digits=2)

#set working directory to folder that contains the data
setwd()


#read in data
data <- read.dta("BES200506080910.dta")

############################################
### Aggregate changes in feeling thermometer
### Figure 1 paper
############################################

labour_ft_pre_2005 <- round(mean(data$labour_ft_pre_2005, na.rm = TRUE), digits = 2)
labour_ft_post_2005 <- round(mean(data$labour_ft_post_2005, na.rm = TRUE),  digits = 2)

cons_ft_pre_2005 <- round(mean(data$cons_ft_pre_2005, na.rm = TRUE),  digits = 2)
cons_ft_post_2005 <- round(mean(data$cons_ft_post_2005, na.rm = TRUE),  digits = 2)

libdem_ft_pre_2005 <- round(mean(data$libdem_ft_pre_2005, na.rm = TRUE),  digits = 2)
libdem_ft_post_2005 <- round(mean(data$libdem_ft_post_2005, na.rm = TRUE),  digits = 2)

labour_ft_pre_2010 <- round(mean(data$labour_ft_pre_2010, na.rm = TRUE),  digits = 2)
labour_ft_post_2010 <- round(mean(data$labour_ft_post_2010, na.rm = TRUE),  digits = 2)

cons_ft_pre_2010 <- round(mean(data$cons_ft_pre_2010, na.rm = TRUE),  digits = 2)
cons_ft_post_2010 <- round(mean(data$cons_ft_post_2010, na.rm = TRUE),  digits = 2)

libdem_ft_pre_2010 <- round(mean(data$libdem_ft_pre_2010, na.rm = TRUE),  digits = 2)
libdem_ft_post_2010 <- round(mean(data$libdem_ft_post_2010, na.rm = TRUE),  digits = 2)

time <- rep(c("pre","post"), 6)
party <- rep(c("labour","labour","cons","cons","libdem","libdem"), 2)
year <- c(rep(2005,6), rep(2010,6))
aggregate <- data.frame(time, party, year)
aggregate$eval <- c(labour_ft_pre_2005, labour_ft_post_2005, cons_ft_pre_2005, cons_ft_post_2005, libdem_ft_pre_2005, libdem_ft_post_2005, labour_ft_pre_2010, labour_ft_post_2010, cons_ft_pre_2010, cons_ft_post_2010, libdem_ft_pre_2010, libdem_ft_post_2010)
aggregate$time <- factor(aggregate$time, levels(aggregate$time)[c(2,1)])

p <- ggplot(aggregate[1:6,], aes(x=time, y=eval, group=party, linetype = party))
p <- p + geom_line() + geom_point() + geom_text(aes(label=eval), hjust=-.20, vjust=0.75)
p <- p + scale_x_discrete("2005 Elections") + scale_y_continuous("Average Feeling Thermometer", limits = c(0, 6)) + theme_bw()
p <- p + theme(axis.title.x = element_text(face="plain", color="black", size=15, hjust = 0.5,vjust = 0.5))
p <- p + theme(axis.title.y = element_text(face="plain", color="black", size=15, hjust = 0.5, vjust = 0.5))
p <- p + theme(legend.title=element_blank()) + theme(legend.text = element_text(size = 15))
print(p)

p <- ggplot(aggregate[7:12,], aes(x=time, y=eval, group=party, linetype = party))
p <- p + geom_line() + geom_point() + geom_text(aes(label=eval), hjust=-.20, vjust=0.75)
p <- p + scale_x_discrete("2010 Elections") + scale_y_continuous("Average Feeling Thermometer", limits = c(0, 6)) + theme_bw()
p <- p + theme(axis.title.x = element_text(face="plain", color="black", size=15, hjust = .5,vjust = 0.5))
p <- p + theme(axis.title.y = element_text(face="plain", color="black", size=15, hjust = .5, vjust = 0.5))
p <- p + theme(legend.title=element_blank()) + theme(legend.text = element_text(size = 15))
print(p)

###########################################
#merge opinion data with media content data
###########################################

#rename variables

data$outlet.05 <- as.character(data$pre_q147)
data$outlet.10 <- as.character(data$aaq147)

media.05 <- read.dta("vistone2005.dta")
media.10 <- read.dta("vistone2010.dta")

#subset data to include only readers of newspaper for which coded content is available

newspapers <- unique(media.05$newspaper)
data <- data[data$pre_q147 %in% newspapers & data$aaq147 %in% newspapers ,]

media.05 <- media.05 %>%
    group_by(newspaper) %>%
    summarise(vislabour  = sum(vislabour, na.rm = TRUE), viscon  = sum(viscon, na.rm = TRUE), vislibdem  = sum(vislibdem, na.rm = TRUE), vislabour_relative = mean(vislabour_relative, na.rm = TRUE), viscon_relative = mean(viscon_relative, na.rm = TRUE), vislibdem_relative= mean(vislibdem_relative, na.rm = TRUE), tonelabour  = sum(tonelabour, na.rm = TRUE), tonecon  = sum(tonecon, na.rm = TRUE), tonelibdem  = sum(tonelibdem, na.rm = TRUE), tonelabour_relative = mean(tonelabour_relative, na.rm = TRUE), tonecon_relative = mean(tonecon_relative, na.rm = TRUE), tonelibdem_relative= mean(tonelibdem_relative, na.rm = TRUE))

media.05 <- as.data.frame(media.05)

media.10 <- media.10 %>%
    group_by(newspaper) %>%
    summarise(vislabour  = sum(vislabour, na.rm = TRUE), viscon  = sum(viscon, na.rm = TRUE), vislibdem  = sum(vislibdem, na.rm = TRUE), vislabour_relative = mean(vislabour_relative, na.rm = TRUE), viscon_relative = mean(viscon_relative, na.rm = TRUE), vislibdem_relative= mean(vislibdem_relative, na.rm = TRUE), tonelabour  = sum(tonelabour, na.rm = TRUE), tonecon  = sum(tonecon, na.rm = TRUE), tonelibdem  = sum(tonelibdem, na.rm = TRUE), tonelabour_relative = mean(tonelabour_relative, na.rm = TRUE), tonecon_relative = mean(tonecon_relative, na.rm = TRUE), tonelibdem_relative= mean(tonelibdem_relative, na.rm = TRUE))

media.10 <- as.data.frame(media.10)


#plug in visibility and tone variables to media content data

for (i in 1:length(data)) {
    if (is.na(data$outlet.05[i]) | is.na(data$outlet.10[i]) ) next
  
    data$vislabour05[i] <- media.05$vislabour[which(media.05$newspaper == data$outlet.05[i])]
    data$vislabour10[i] <- media.10$vislabour[which(media.10$newspaper == data$outlet.10[i])]
    data$viscons05[i] <- media.05$viscon[which(media.05$newspaper == data$outlet.05[i])]
    data$viscons10[i] <-media.10$viscon[which(media.10$newspaper == data$outlet.10[i])]
    data$vislibdem05[i] <- media.05$vislibdem[which(media.05$newspaper == data$outlet.05[i])]
    data$vislibdem10[i] <-media.10$vislibdem[which(media.10$newspaper == data$outlet.10[i])]

    data$tonelabour05[i] <- media.05$tonelabour[which(media.05$newspaper == data$outlet.05[i])]
    data$tonelabour10[i] <-media.10$tonelabour[which(media.10$newspaper == data$outlet.10[i])]
    data$tonecons05[i] <- media.05$tonecon[which(media.05$newspaper == data$outlet.05[i])]
    data$tonecons10[i] <-media.10$tonecon[which(media.10$newspaper == data$outlet.10[i])]
    data$tonelibdem05[i] <- media.05$tonelibdem[which(media.05$newspaper == data$outlet.05[i])]
    data$tonelibdem10[i] <-media.10$tonelibdem[which(media.10$newspaper == data$outlet.10[i])]

}

#drop variables
data$outlet.05 <- data$outlet.10 <- NULL

#recode feeling thermometers so that positive scores denote positive change
data$evalslabour05 <- -1*data$labour_ft_delta_05
data$evalslabour10 <- -1*data$labour_ft_delta_10
data$evalscons05 <- -1*data$cons_ft_delta_05
data$evalscons10 <- -1*data$cons_ft_delta_10
data$evalslibdem05 <- -1*data$libdem_ft_delta_05
data$evalslibdem10 <- -1*data$libdem_ft_delta_05

#recode visibility and tone variables
data$zvislabour05 <- data$vislabour05 /20
data$zviscons05 <- data$viscons05 /20
data$zvislibdem05 <- data$vislibdem05 /20
data$ztonelabour05 <- data$tonelabour05 /20
data$ztonecons05 <- data$tonecons05 /20
data$ztonelibdem05 <- data$tonelibdem05 /20

data$zvislabour10 <- data$vislabour10 /20
data$zviscons10 <- data$viscons10 /20
data$zvislibdem10 <- data$vislibdem10 /20
data$ztonelabour10 <- data$tonelabour10 /20
data$ztonecons10 <- data$tonecons10 /20
data$ztonelibdem10 <- data$tonelibdem10 /20


########################################
#Figure Descriptive Statistics: Appendix
########################################

newspaper.data <- data
newspaper.data <- newspaper.data[order(newspaper.data$pre_q147),]

vis.tone.05 <- ddply(newspaper.data, .(pre_q147), summarise, labourvis = round(mean(vislabour05, na.rm = TRUE), 0), labourtone = round(mean(tonelabour05, na.rm = TRUE), 0),consvis = round(mean(viscons05, na.rm = TRUE), 0), constone = round(mean(tonecons05, na.rm = TRUE), 0), libdemvis = round(mean(vislibdem05, na.rm = TRUE), 0), libdemtone = round(mean(tonelibdem05, na.rm = TRUE), 0))

names(vis.tone.05) <- c("newspaper", "lab vis", "lab tone", "cons vis", "cons tone", "lib dem vis", "lib dem tone")
vis.tone.05 <- vis.tone.05[-c(2, 9), ]
rownames(vis.tone.05) <- NULL

vis.tone.10 <- ddply(newspaper.data, .(pre_q147), summarise, labourvis = round(mean(vislabour10, na.rm = TRUE), 0), labourtone = round(mean(tonelabour10, na.rm = TRUE), 0),consvis = round(mean(viscons10, na.rm = TRUE), 0), constone = round(mean(tonecons10, na.rm = TRUE), 0), libdemvis = round(mean(vislibdem10, na.rm = TRUE), 0), libdemtone = round(mean(tonelibdem10, na.rm = TRUE), 0))

names(vis.tone.10) <- c("newspaper", "lab vis", "lab tone", "cons vis", "cons tone", "lib dem vis", "lib dem tone")
vis.tone.10 <- vis.tone.10[-c(2, 9), ]
rownames(vis.tone.10) <- NULL

output.05 <- stargazer(vis.tone.05, 
                    type = "html", no.space=TRUE, single.row = F,
                    summary=FALSE,  model.names = FALSE, dep.var.labels.include = FALSE, dep.var.caption = "XXX",
                    notes.append = FALSE)

output.10 <- stargazer(vis.tone.10, 
                    type = "html", no.space=TRUE, single.row = F,
                    summary=FALSE,
                    notes.append = FALSE)

############################################
#Figure Descriptive Statistics: not in paper
############################################

newspaper.data <- subset(data, selected == 1 & selected10 == 1)
newspaper.data <- subset(newspaper.data, as.numeric(newspaper.data$pre_q146) == 1 | as.numeric(newspaper.data$pre_q146) == 2 | as.numeric(newspaper.data$aaq145) == 1 | as.numeric(newspaper.data$aaq145) == 2)

variables <- c("lowed", "meded", "hied", "gender", "income", "partyid_labour", "partyid_cons", "partyid_libdem", "white")

newspaper.data <- newspaper.data[,variables]

names(newspaper.data) <- c("Low Education", "Medium Education", "High Education", "Gender", "Income", "Party Identification", "Conservative Identification", "Lib Dem Identification", "Race")

stargazer(newspaper.data[,variables], digits = 2)

output <- stargazer(newspaper.data, 
                    type = "html", title = "Summary Statistics", no.space=TRUE, single.row = F,
                    notes.append = FALSE, digits = 2)

#read in function to create correlation table
source("corstarsl.R")

output <- print(xtable(corstarsl(newspaper.data)), type = "html")


#########################################
#Within-Elections, between-subjects model
#Table 2 paper
#########################################

#select variables

variables <- c("besid", "evalslabour05", "evalscons05", "evalslibdem05", "evalslabour10",  "evalscons10", "evalslibdem10", "zvislabour05", "zviscons05", "zvislibdem05","zvislabour10", "zviscons10", "zvislibdem10",  "ztonelabour05", "ztonecons05", "ztonelibdem05","ztonelabour10", "ztonecons10", "ztonelibdem10", "partyid_labour", "partyid_cons", "partyid_libdem", "lowed", "meded", "gender", "income", "white")

model.data <- data[,variables]
model.data <- na.omit(model.data)

data.long <- model.data
data.long <- data.long[,variables]

names(data.long) <- c("id", "evals.1", "evals.2", "evals.3", "evals.4", "evals.5", "evals.6", "vis.1", "vis.2", "vis.3", "vis.4", "vis.5", "vis.6", "tone.1", "tone.2", "tone.3", "tone.4", "tone.5", "tone.6","partyid_labour", "partyid_cons", "partyid_libdem", "lowed", "meded", "gender", "income", "white")

data.long <-reshape(data.long, varying= c(evals = names(data.long)[2:7], visibility = names(data.long)[8:13], tone = names(data.long[14:19])), direction= "long", idvar= names(data.long)[1])

coefficients.model <- c("Visibility", "Tone", "Labour Identification", "Conservative Identification", "Liberal Democrat Identification", "Low Education", "Medium Education", "Gender", "Income", "Race")


#estimate within-elections, between-subjects model
model <- lm(evals ~ vis + tone + partyid_labour + partyid_cons + partyid_libdem + lowed + meded + gender + income + white, data = data.long)
zelig.between.subjects.within.elections <- zelig(evals ~ vis + tone + partyid_labour + partyid_cons + partyid_libdem + lowed + meded + gender + income + white, model = "normal", data = data.long)

output <- stargazer(model, 
                    type = "html", title = "Between-Subjects, Within-Elections Comparisons 2005", no.space=TRUE, single.row = T,
                    covariate.labels = coefficients.model, initial.zero = T,
                    model.names = FALSE, dep.var.labels.include = FALSE, dep.var.caption = "XXX",
                    notes = "Standardized coefficients, standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.10",
                    style = "io",
                    star.cutoffs = c(0.10, 0.05, 0.01),
                    star.char = c("*", "**", "***"),
                    notes.append = FALSE)

output <- output[-c(4)]


##############################################################
#Between-subjects, within-elections treatment effects figure 2
##############################################################

N <- 1
first.differences <- data.frame(
    party = character(N),
    mean.change.vis = numeric(N),
    sd.change.vis = numeric(N),
    mean.change.tone = numeric(N),
    sd.change.tone = numeric(N),
    stringsAsFactors = FALSE)

first.differences$party <- c("Party")


#######################
# Visibility
#######################

x.low <- setx(zelig.between.subjects.within.elections, vis = min(data.long$vis), fn = list(numeric = mfv))
x.high <- setx(zelig.between.subjects.within.elections, vis = max(data.long$vis), fn = list(numeric = mfv))

s.out.zelig.between.subjects.within.elections <- sim(zelig.between.subjects.within.elections, x = x.low, x1 = x.high, num = 100000) 

qi <- zelig_qi_to_df(s.out.zelig.between.subjects.within.elections)

first.differences[1,2] <- mean(qi$expected_value[qi$setx_value == "x1" ]) - mean(qi$expected_value[qi$setx_value == "x" ])
first.differences[1,3] <- sd(qi$expected_value[qi$setx_value == "x1" ] - qi$expected_value[qi$setx_value == "x" ])


#######################
# tone
#######################

x.low <- setx(zelig.between.subjects.within.elections, tone = min(data.long$tone), fn = list(numeric = mfv))
x.high <- setx(zelig.between.subjects.within.elections, tone = max(data.long$tone), fn = list(numeric = mfv))

s.out.zelig.between.subjects.within.elections <- sim(zelig.between.subjects.within.elections, x = x.low, x1 = x.high, num = 100000) 

qi <- zelig_qi_to_df(s.out.zelig.between.subjects.within.elections)

first.differences[1,4] <- mean(qi$expected_value[qi$setx_value == "x1" ]) - mean(qi$expected_value[qi$setx_value == "x" ])
first.differences[1,5] <- sd(qi$expected_value[qi$setx_value == "x1" ] - qi$expected_value[qi$setx_value == "x" ])

#generate visibility plot
vis <- ggplot(first.differences,
              aes(x = party , y = mean.change.vis,
                  ymin = mean.change.vis - 2*sd.change.vis,
                  ymax = mean.change.vis + 2*sd.change.vis))
vis <- vis + geom_pointrange() + geom_abline(intercept = 0, slope = 0, color = "black") + expand_limits(y=0)
vis <- vis + theme_minimal() + scale_y_continuous("Effect Size") + theme(plot.title = element_text(hjust = 0.5))
vis <- vis + theme(axis.title.y = element_text(face="bold", color="black", size=15, hjust = .5, vjust = 1)) + theme(axis.title.x = element_blank())
vis <- vis + ggtitle("First Differences, Visiblity") + theme(plot.title = element_text(lineheight=.8, face="bold"))
print(vis)

#generate tone plot
tone <- ggplot(first.differences,
              aes(x = party , y = mean.change.tone,
                  ymin = mean.change.tone - 2*sd.change.tone,
                  ymax = mean.change.tone + 2*sd.change.tone))
tone <- tone + geom_pointrange() + geom_abline(intercept = 0, slope = 0, color = "black") + expand_limits(y=0)
tone <- tone + theme_minimal() + scale_y_continuous("Effect Size") + theme(plot.title = element_text(hjust = 0.5))
tone <- tone + theme(axis.title.y = element_text(face="bold", color="black", size=15, hjust = .5, vjust = 1)) + theme(axis.title.x = element_blank())
tone <- tone + ggtitle("First Differences, Tone") + theme(plot.title = element_text(lineheight=.8, face="bold"))
print(tone)

###########################################
#Within-subjects, within-elections analysis
#Table 3 in the paper
###########################################

data.long <- model.data

variables <- c("besid", "evalslabour05", "evalscons05", "evalslibdem05", "evalslabour10",  "evalscons10", "evalslibdem10", "zvislabour05", "zviscons05", "zvislibdem05","zvislabour10", "zviscons10", "zvislibdem10",  "ztonelabour05", "ztonecons05", "ztonelibdem05","ztonelabour10", "ztonecons10", "ztonelibdem10")

data.long <- data.long[,variables]

names(data.long) <- c("id", "evals.1", "evals.2", "evals.3", "evals.4", "evals.5", "evals.6", "vis.1", "vis.2", "vis.3", "vis.4", "vis.5", "vis.6", "tone.1", "tone.2", "tone.3", "tone.4", "tone.5", "tone.6")

data.long <-reshape(data.long, varying= c(evals = names(data.long)[2:7], visibility = names(data.long)[8:13], tone = names(data.long[14:19])), direction= "long", idvar= names(data.long)[1])
                   
coefficients <- c("Visibility", "Tone")

#estimate within-subjects, within-elections model
within.subjects.within.elections <- lm(evals ~ vis + tone + factor(id), data=data.long)
zelig.within.subjects.within.elections <- zelig(evals ~ vis + tone + factor(id), model = "normal", data = data.long)

output <- stargazer(within.subjects.within.elections, 
                    type = "html", no.space=TRUE, single.row = T,
                    column.labels = "Change in Party Evaluations",
                    covariate.labels = coefficients,
                    model.names = FALSE, dep.var.labels.include = FALSE, dep.var.caption = "XXX",
                    notes = "Standardized coefficients, standard errors in parentheses. *** p<0.001, ** p<0.01, * p<0.05",
                    style = "io",
                    star.cutoffs = c(0.05, 0.01, 0.001),
                    star.char = c("*", "**", "***"),
                    notes.append = FALSE)

output <- output[c(1:4, 557:563)]


#############################################################
#Within-subjects, within-elections treatment effects figure 3
#############################################################

N <- 1
first.differences <- data.frame(
    party = character(N),
    mean.change.vis = numeric(N),
    sd.change.vis = numeric(N),
    mean.change.tone = numeric(N),
    sd.change.tone = numeric(N),
    stringsAsFactors = FALSE)

first.differences$party <- c("Party")


#######################
# Visibility
#######################

x.low <- setx(zelig.within.subjects.within.elections, vis = min(data.long$vis), fn = list(numeric = mfv))
x.high <- setx(zelig.within.subjects.within.elections, vis = max(data.long$vis), fn = list(numeric = mfv))

s.out.zelig.within.subjects.within.elections <- sim(zelig.within.subjects.within.elections, x = x.low, x1 = x.high, num = 100000) 

qi <- zelig_qi_to_df(s.out.zelig.within.subjects.within.elections)

first.differences[1,2] <- mean(qi$expected_value[qi$setx_value == "x1" ]) - mean(qi$expected_value[qi$setx_value == "x" ])
first.differences[1,3] <- sd(qi$expected_value[qi$setx_value == "x1" ] - qi$expected_value[qi$setx_value == "x" ])


#######################
# tone
#######################

x.low <- setx(zelig.within.subjects.within.elections, tone = min(data.long$tone), fn = list(numeric = mfv))
x.high <- setx(zelig.within.subjects.within.elections, tone = max(data.long$tone), fn = list(numeric = mfv))

s.out.zelig.within.subjects.within.elections <- sim(zelig.within.subjects.within.elections, x = x.low, x1 = x.high, num = 100000) 

qi <- zelig_qi_to_df(s.out.zelig.within.subjects.within.elections)

first.differences[1,4] <- mean(qi$expected_value[qi$setx_value == "x1" ]) - mean(qi$expected_value[qi$setx_value == "x" ])
first.differences[1,5] <- sd(qi$expected_value[qi$setx_value == "x1" ] - qi$expected_value[qi$setx_value == "x" ])

#generate visibility plot
vis <- ggplot(first.differences,
              aes(x = party , y = mean.change.vis,
                  ymin = mean.change.vis - 2*sd.change.vis,
                  ymax = mean.change.vis + 2*sd.change.vis))
vis <- vis + geom_pointrange() + geom_abline(intercept = 0, slope = 0, color = "black") + expand_limits(y=0)
vis <- vis + theme_minimal() + scale_y_continuous("Effect Size") + theme(plot.title = element_text(hjust = 0.5))
vis <- vis + theme(axis.title.y = element_text(face="bold", color="black", size=15, hjust = .5, vjust = 1)) + theme(axis.title.x = element_blank())
vis <- vis + ggtitle("First Differences, Visiblity") + theme(plot.title = element_text(lineheight=.8, face="bold"))
print(vis)

#generate tone plot
tone <- ggplot(first.differences,
              aes(x = party , y = mean.change.tone,
                  ymin = mean.change.tone - 2*sd.change.tone,
                  ymax = mean.change.tone + 2*sd.change.tone))
tone <- tone + geom_pointrange() + geom_abline(intercept = 0, slope = 0, color = "black") + expand_limits(y=0)
tone <- tone + theme_minimal() + scale_y_continuous("Effect Size") + theme(plot.title = element_text(hjust = 0.5))
tone <- tone + theme(axis.title.y = element_text(face="bold", color="black", size=15, hjust = .5, vjust = 1)) + theme(axis.title.x = element_blank())
tone <- tone + ggtitle("First Differences, Tone") + theme(plot.title = element_text(lineheight=.8, face="bold"))
print(tone)

####################################
#Within-subjects, between-elections
#Table 4 in the paper
####################################

data.long <- model.data
data.long$d.evals.labour <- data.long$evalslabour10 - data.long$evalslabour05
data.long$d.evals.cons <- data.long$evalscons10 - data.long$evalscons05
data.long$d.evals.libdem <- data.long$evalslibdem10 - data.long$evalslibdem05

data.long$d.zvis.labour <- data.long$zvislabour10 - data.long$zvislabour05
data.long$d.zvis.cons <- data.long$zviscons10 - data.long$zviscons05
data.long$d.zvis.libdem <- data.long$zvislibdem10 - data.long$zvislibdem05

data.long$d.ztone.labour <- data.long$ztonelabour10 - data.long$ztonelabour05
data.long$d.ztone.cons <- data.long$ztonecons10 - data.long$ztonecons05
data.long$d.ztone.libdem <- data.long$ztonelibdem10 - data.long$ztonelibdem05

variables <- c("besid", "d.evals.labour", "d.evals.cons", "d.evals.libdem", "d.zvis.labour", "d.zvis.cons", "d.zvis.libdem", "d.ztone.labour", "d.ztone.cons", "d.ztone.libdem")

data.long <- data.long[,variables]

names(data.long) <- c("id", "evals.1", "evals.2", "evals.3", "vis.1", "vis.2", "vis.3", "tone.1", "tone.2", "tone.3")

data.long <-reshape(data.long, varying= c(evals = names(data.long)[2:4], visibility = names(data.long)[5:7], tone = names(data.long[8:10])), direction= "long", idvar= names(data.long)[1])


#estimate within-subjects, between elections model
within.subjects.between.elections <- lm(as.numeric(evals) ~ vis + tone + factor(id) - 1, data= data.long)
zelig.within.subjects.between.elections <- zelig(evals ~ vis + tone + factor(id), model = "normal", data = data.long)


coefficients <- c("Visibility", "Tone")

output <- stargazer(within.subjects.between.elections, 
                    type = "html", no.space=TRUE, single.row = T,
                    column.labels = "Difference in Change in Party Evaluations",
                    covariate.labels = coefficients,
                    model.names = FALSE, dep.var.labels.include = FALSE, dep.var.caption = "XXX",
                    notes = "Standardized coefficients, standard errors in parentheses. *** p<0.001, ** p<0.01, * p<0.05",
                    style = "io",
                    star.cutoffs = c(0.05, 0.01, 0.001),
                    star.char = c("*", "**", "***"),
                    notes.append = FALSE)

output <- output[c(1:4, 557:563)]

###################################################
#within-subjects, between-elections, figure 4 paper
###################################################

N <- 1
first.differences <- data.frame(
    party = character(N),
    mean.change.vis = numeric(N),
    sd.change.vis = numeric(N),
    mean.change.tone = numeric(N),
    sd.change.tone = numeric(N),
    stringsAsFactors = FALSE)

first.differences$party <- c("Party")


#######################
# Visibility
#######################

x.low <- setx(zelig.within.subjects.between.elections, vis = min(data.long$vis), fn = list(numeric = mfv))
x.high <- setx(zelig.within.subjects.between.elections, vis = max(data.long$vis), fn = list(numeric = mfv))

s.out.zelig.within.subjects.between.elections <- sim(zelig.within.subjects.between.elections, x = x.low, x1 = x.high, num = 100000) 

qi <- zelig_qi_to_df(s.out.zelig.within.subjects.between.elections)

first.differences[1,2] <- mean(qi$expected_value[qi$setx_value == "x1" ]) - mean(qi$expected_value[qi$setx_value == "x" ])
first.differences[1,3] <- sd(qi$expected_value[qi$setx_value == "x1" ] - qi$expected_value[qi$setx_value == "x" ])


#######################
# tone
#######################
x.low <- setx(zelig.within.subjects.between.elections, tone = min(data.long$tone), fn = list(numeric = mfv))
x.high <- setx(zelig.within.subjects.between.elections, tone = max(data.long$tone), fn = list(numeric = mfv))

s.out.zelig.within.subjects.between.elections <- sim(zelig.within.subjects.between.elections, x = x.low, x1 = x.high, num = 100000) 

qi <- zelig_qi_to_df(s.out.zelig.within.subjects.between.elections)

first.differences[1,4] <- mean(qi$expected_value[qi$setx_value == "x1" ]) - mean(qi$expected_value[qi$setx_value == "x" ])
first.differences[1,5] <- sd(qi$expected_value[qi$setx_value == "x1" ] - qi$expected_value[qi$setx_value == "x" ])


#generate visibility plot

vis <- ggplot(first.differences,
              aes(x = party , y = mean.change.vis,
                  ymin = mean.change.vis - 2*sd.change.vis,
                  ymax = mean.change.vis + 2*sd.change.vis))
vis <- vis + geom_pointrange() + geom_abline(intercept = 0, slope = 0, color = "black") + expand_limits(y=0)
vis <- vis + theme_minimal() + scale_y_continuous("Effect Size") + theme(plot.title = element_text(hjust = 0.5))
vis <- vis + theme(axis.title.y = element_text(face="bold", color="black", size=15, hjust = .5, vjust = 1)) + theme(axis.title.x = element_blank())
vis <- vis + ggtitle("First Differences, Visiblity") + theme(plot.title = element_text(lineheight=.8, face="bold"))
print(vis)

#generate tone plot

tone <- ggplot(first.differences,
              aes(x = party , y = mean.change.tone,
                  ymin = mean.change.tone - 2*sd.change.tone,
                  ymax = mean.change.tone + 2*sd.change.tone))
tone <- tone + geom_pointrange() + geom_abline(intercept = 0, slope = 0, color = "black") + expand_limits(y=0)
tone <- tone + theme_minimal() + scale_y_continuous("Effect Size") + theme(plot.title = element_text(hjust = 0.5))
tone <- tone + theme(axis.title.y = element_text(face="bold", color="black", size=15, hjust = .5, vjust = 1)) + theme(axis.title.x = element_blank())
tone <- tone + ggtitle("First Differences, Tone") + theme(plot.title = element_text(lineheight=.8, face="bold"))
print(tone)

#END
